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ABSTRACT 


Using linear systems theory as a framework, the solution for the acoustic field 
present in a range-independent acoustic channel excited by a complex -weighted, 
planar array of point sources with an arbitrary input electrical signal is derived. 
The ocean medium is characterized by a transfer function, obtainable as the solution 
to the Helmholtz wave equation. The transfer function for an isospeed, three-layer 
waveguide is derived. The unbounded homogeneous medium equations are derived 
as a special case of the waveguide problem. The problem of interference due to the 
presence of a pressure-release surface is also derived as a special case. The linear 
systems approach lends itself to a modular computer implementation, in which 
different ocean medium models are represented by subroutines implementing their 
transfer functions. The equations for a range-independent medium are implemented 
as a group of subprograms. Results are presented for the special cases of a 
homogeneous medium and the surface reflection problem, which can be checked 
against known, easily interpreted analytical solutions. Finally, an example of 


waveform prediction for the isospeed, three-layer waveguide is presented. 
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I. INTRODUCTION 


Model-based or matched-field signal processing refers to the use of the 
knowledge of the physics of a problem for the construction of mathematical models 
and its application to suitable signal processing algorithms. Its application to 
underwater acoustics has been the subject of many papers in the signal processing 
literature, as in Ziomek and Blount [Ref. 1], Baggeroer, Kuperman and 
Schmidt [Ref. 2] and references therein. As the interest in this field grows and new 
algorithms are developed, so will the need for simulation tools to verify the behavior 
of those algorithms. Such tools should be able to generate signals with spatial and 
temporal structures like those found in real environments. 

The waveform prediction problem has been studied by Officer [Ref. 3:pp. 
101-110,130] for the specific cases of pulse reflection from a boundary and 
transmission in shallow water. DiNapoli and Deavenport [Ref.4:pp. 131-134] 
describe briefly a method employing the Fast-Field-Program (FFP) suitable to 
range-independent problems. Our purpose is to derive a generic solution for the 
range-independent. time-invariant, deterministic acoustic channel excited by an 
array of point sources transmitting an arbitrary signal. In our derivation we will 
use a linear systems approach to the acoustic problem, as described in 
Ziomek [Ref. 5], whose notation we follow closely. In particular, we will apply our 
results to the isospeed Pekeris waveguide. As we will show, the linear systems 
approach lends itself to a modular computer implementation, where different ocean 
media are represented by subroutines or functions implementing their transfer 


function. 


In Section Il we present the basic input-output relationships for a 
space-variant, time-independent linear system, which are the basis of our 
derivations. Next, we show how those equations can be used to represent an 
acoustic channel. The particular case of a range-independent ocean is then studied, 
when we derive the output waveform equation and describe the process of derivation 
of the medium transfer function. Finally, we apply the above results to the isospeed 
Pekeris waveguide. In Section III we discuss computer implementation issues and 
present some results of a computer implementation of the Pekeris waveguide 


equations. 


II. ANALYSIS 


A. THE ACOUSTIC CHANNEL AS A LINEAR SYSTEM 
In this section we will present some results from the theory of space-time linear 
systems [Ref. 5], establish the notation to be used, and derive some intermediate 
results to be used in subsequent sections. 
1. Space-Variant, Time-Invariant Linear Systems 
a. Impulse Response 
Linear, time-invariant, space-variant filters, as shown in Fig. 1, are 


represented by linear partial differential equations whose coefficients are time 
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independent. In linear systems terminology, suc! 
time-invariant, space-dependent impulse response h(7,ro; r). The output y(tr) is 


given by [Ref. 5;p. 3] 


(tr) = Í Í Ol ar I, i ial andre. (1) 


00 00 
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Fig. 1. A Linear, Time-Invariant, Space-Variant System 


A physical interpretation of A, 7 and r, can be obtained by letting the 
input be an impulse applied at time t= and at position r = Ts, that is, 
g(tr) = ó(t-tor-rs) — é(t-t5) é(r-rs). Substituting this expression for g(t,r) into 
Eq. (1) and using the sifting property of the impulse function we obtain, for the 


output, 


EET T (2) 
v 
T [a 


from which we can define h(r,ro;r) as the response of the system at time t and 
position r due to the application of an impulse 7 seconds ago, that is, at the instant 
tb = t-7, when the point source is positioned at rs = r-ro. Figure 2 illustrates the 
geometry of this problem. 
b. Transfer Function 
Let the input to the system be a plane wave, that is, g(tr) = 
- drv tr) with frequency $ and spatial frequency vector vo. The output of 


the system is given by Eq. (1) as 


ptr) - f E ro) hU oro r (3a) 


p( tr) = on | f TARN oT) a (8h) 
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Basic Geometry for the Interpretation of hl r.r.:r). 


OI 


(tr) e d? rg. A ej (3c) 


T- 
V= o 


where 8. pTepresents the time-domain Fourier transform and Sr _, Fepresents the 
E O 

spatial-domain Fourier transform. The subscripts under the symbol $ represent the 

variables involved in each transform. We define H(fwv;r), the system transfer 


function, as follows: 
Ház R 1 m "dp -J 2v( fr-v-ro) E 
A dE uA hr roi) j = e h(T,ro;r) dr dr; . (4) 


Using this definition, the output of the system can be written as 
d - 
p( hr) = H( froin) ef “Mot vot) (5) 


the expected result from linear systems theory. 

The parameter ro in the integrations in Eqs. (1) and (4) represents the 
vectorial difference between the observation point r and the point source position rs. 
The integrations must be performed taking r as a constant, that is, by changing the 
source position rs so that ro assumes all possible values in R3, and the results, y(ér) 


or H(f,vo;r), are valid for that ''fixed" observation point r. 
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c. Response to a Time-Harmonic Point Source 
When evaluating the transfer function, it is convenient to use a 
time -harmonic point source as the input, g(tr) = ei Tht ó(r-rs) with frequency $ 


and location rs. The output, from Eqs. (1) and (4), is given by 


p(ir) = Í — ó(r-rg-ro) h(m,ro;r) dr dro, (6a) 
CEU Jed PLE h(T,r-rg;r) dr, (6b) 





alt) = his (bl) = po (6c) 
Io= I-Is 
OI 
(ina dira! {Ho}. (6d) 


d. Transform Relationships 
(1) Fourier Transform Definitions and Notation. The time-domain 
Fourier transform of a function g(tr) = g(t,x,y,z) and its inverse are given. 


respectively, by 


GUI) = By Lot Yo [| an ena 7 


=] 


and 


g(6r) = #1 pf Gin) } = Í GC) DIN lap. (8) 


00 


The spatial Fourier transform and its inverse are given, respectively, by 


G(tv) = 8,_, { 940) } =f a(t) AT de, (9) 
and 
g(tr) =F! { (tv) } = {ct IR Ty | (10) 


o 


where dr = dr dy dz , dv = d& d dk . and wvr=tk+yh+z£. The above 
transforms involve triple integrals and can be split into three spatial transforms in 7 


(£), v (f) and z (£). For example, the transforms in 7 (f) are written as 


~ i no 7 
Gb bow) = 8, ¢ { (6542) } = | (tx y) e) "tds, (11) 


To 


and 


(n9) — S, e CGU o2) - Í Cli kna) e a (2) 


oo 


The transforms in y (f) and z (£) follow the same notation. The full space-time 


transforms can be written as 


Cf) =p AD) =p [Ep ES, eS p Co), (3a) 


or 
G(fv) = rea g(tr) dt dr, (13b) 
and 
gtr) = 8). RL { GY) = pp Bop (OU), (Ida) 
or 


g(t.r) - [ y JJ 2 T tov) G(fv) df dv . (14b) 


The above relationships are valid for the impulse response/transfer function pair if ! 


is replaced by 7r and r is replaced by ro = (25, Y,%0), that is, 
A frost) 2 8, pt hino }, (15) 


H(trr) 8 AMA), (16) 


H(t fyo 20;1) = 8, Nx h (7,10, BD) }, (17) 
and 
GA o h(m,ro;r) ) . (4) 


The same applies for the respective inverse transforms. 
(2) Input-Output Relationships. During our developments, we will 
need to express the output of the system, y(t,r), in terms of the transformed input 


function G( fy) and the transfer function. From the Fourier transform definition 
g(tr) = f Í PANETE (18) 


follows the transform property 


g(t-7,r-1.) = Í Í Jemft-ver) II df dv. (19) 


Substituting this expression into Eg. (1) yields 


p(tr) = Bee Pe Arr vote) gp gy) 


x h(m,ro;r) df dv dr dro, (20a) 
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p(tr) = Í [urna DG in| Jer SMT o qeeedy, (20b) 


HL v:r) 


" f J “G( fu) H( frre) d 27 mapa. (20c) 


which is the desired relationship. 

Analogous to a time-variant system, for which the frequencies at 
the output are different from the input, a space-variant system causes the output 
spatial frequency vector to be different from the input. We will use 8 to represent 
the output spatial frequency vector. The input-output relationship in the 
transformed form can be obtained by substituting Eq. (20c) into the Fourier 


transform definition 
jun | | e a de (21) 


Doing so vields 


}(fB) = ff] IE H(n,v;r) eni give. t) dv DÀ 2r CB qu dr, (22a) 


(f) = [f feno H(n,v;r) 23 27(v-B)-1 f 


00 


IE a dr dv, (29b) 
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(ff) = Í J {f G(nv) Hir) e 3 2 07P s pg in} dedo (22c) 


b( fp) = Í J G(£») H(fy;r) e) 29 9-D ae dy. (22d) 


2. The Acoustic Channel 
As long as we deal with small-amplitude acoustic signals, the governing 
wave equation is linear. We will restrict ourselves to time-invariant problems, 
where the properties of the medium do not change with time, and the source, 
receiver and the boundaries are motionless. 
Figure 3 is a pictorial representation of the acoustic channel. A source or 


" an electrical signal to the medium and a 


transmitter rectangular array "couples 
receiver array transforms the acoustic signal back to electrical form. Both arrays 
are assumed to act as linear filters. We will not be concerned here with the 
generation of the electrical transmitted signal, assumed given, nor with the 
processing done on the received electrical signals at each hydrophone position. The 
acoustic channel can be represented by the block diagram shown in Fig. 1, which 
can be broken down into the three basic components sketched in Fig. 3: transmitter 
array, medium and receiver array. Figure 4 shows these components in block 
diagram form. [Ref. 5;pp. 3-4] 
a. The medium 


The propagation of the acoustic signal through the medium 1s governed 


by the wave equation 


D 


D 
Tz 
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Fig. 3. Geometry of the Acoustic Channel Problem. 
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Fig. 4. Block -Diagram Representation of the Acoustic Channel. 


] m 
v^ e(t ir) - EOD) p(br) ni g(t.r) (23) 
HETI 


and appropriate boundary conditions, where y(tr) is the velocity potential 
aA ) g(tr) is the source term, which represents the rate of mass injection into 
the space per unit mass or rate at which fluid volume is added per unit volume 

ly [Ref. 5;p. 1). These two terms are represented in Fig. 4 as the output and 


input to the acoustic medium block, respectively. The solution of Eq. (23) is 


obtained from Eqs. (1) or (20c) as 


i (cae = f { g( t- T,r-ro) h CTS To;T) disce (24) 


OI 


me - | Ju Eao A (25) 
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The transform of the output (ff) is a function of the spatial 
frequency vector É, which is, in general, different frorn v, due to the space-variant 
nature of the medium. The physical meaning, in terms of ray acoustics, is that a 
non-homogeneous medium (H is a function of r) causes change in the direction of 
propagation of the sound. 

The input -output relationship in the transformed form can be obtained 


from Eq. (22d), resulting in 
(58) = Í | G(r) H (frir) eI TOPE qp ay (26) 


If the system is space-invariant, an unbounded homogeneous medium 
for example, then i is independent of r and Eq. (26) reduces to the expected 


expression for a time-invariant. space-invariant linear system: 


e [ ouo ua | imm ae dy (27a) 
wu e [ou Hu fom àv. (27b) 
DUA = CUA HD (27) 
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b. The Transmitter 
The transmitter array, the first block in Fig. 4, is characterized by a 
complex aperture function An (41) (Amp! sa relating the electrical input z(&r) 


(Amp) to the source distribution g(t,r) (sia through the expression [Ref. 5:p. 30] 


G( fr) = X(fr) A, (Er). (28) 


The far-field directivity function D. (fv) m Ap 370) is related to the complex 


aperture function through the transform [Ref. 5:p. 38] 

D (fv) e 8, LA UD T. (29) 
Using Eq. (29), the spatial Fourier transform of Eq. (28) 1s given by 

GU DIESEN y, DV) (30a) 


OT 


CATE | X( fa) D. (fv-a) da, (30b) 


no 


where the symbol (5) denotes a three-dimensional convolution with respect to ». 
Equations (30a) and (30b) give the relationship between the electrical input and the 
source of the acoustic signal. The relationship between the input electrical signal 


and the velocity potential (the output of the medium block of Fig. 4) is obtained by 
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substituting Eq. (30b) into Eqs. (25) or (26). Let us introduce the simplifying 
assumption that the electrical input signal is independent of position, that is, 
t(tr)= s(t) and X(fv) = X(f) 6(v). This is justified by the fact that the 
complex weights applied to the electrical signal (array shading and beam steering) 
are taken into account in D. (5v), the far-field directivity function. With that 


assumption, Eq. (30a) becomes 
G(4») & XU) $0) * D, (£v), (31a) 
G( fu) = X(f) Dy 40), (31b) 


which, when substituted into Eq.:(25), gives the required relationship between the 


input electrical signal and the velocity potential at the output of the medium block, 
S AM | ME | 
ptr) af | X(f)D (bY) A, (£50) JT vt) dy df. (32) 


A further restriction we impose to our problem is that the transmitter array is 
composed of point sources, has rectangular shape, is parallel to the XY plane, as 
shown in Fig. 3, and is centered at position rs — (%,Ys,%). The array has 


dimensions Me x Ne (first dimension along the X direction and second dimension 


along Y), with interelement spacing dy; meters along X and dy, meters along the Y 
direction. Then, the far-field directivity function for Mẹ and M odd is given 


by (Ref. 5;ch. 4] 


M,-1 N, -1 
Es E 3 


D (v) = 2 >, palf) BE pais (39) 
a M 
ER EE Nido 


where ry4— (pdxt + Ts, qdyt + Ys, 25) represents the position of each element and 
aun = apa Cf ) e pa) is the complex weight associated with the element at 
position [pq 

Cc. The Receiver 
The receiver array, the last block in Fig. 4, is characterized by a 
“Oy 
} 


complex aperture function A, (sr) (Vsm ~} relating the velocity potential g(t.r) 


to the electrical signal y(t,r) (V RA through the expression [Ref. 5:p. 31] 
Vr) e (4) AUN). (34 


As stated earlier, we are not concerned with the signal processing done 
on the received signal. Our purpose is to predict the waveform at the output of 
each receiver array element. So, we consider the term Ap (HI) as pertaining to a 
single transducer element. The quantity y(t,r) is the electrical signal distribution 


(V m due to the velocity potential g(t,r) at a point ron the transducer. In order 
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to obtain the total electrical signal y(t) at the transducer terminals, we must 


integrate y(t,r) over the whole transducer volume, that is, integrate with respect to 


I, 
vif) = J yun dr, 
D= f Un AUD dr (35) 
For a point sensor located at rr, the complex aperture function is given 
by 


A, (fr) = A(f) 6(r-rr) , (36) 


which results in an electrical signal at the terminals given [ see Eq. (35) ] by 


Y(f) 2 €(fr) A(f). (37) 


If the transducer is a pressure sensor, then A( f) is proportional to the frequency f, 
that is, A(f)~jf and Y(f)~ jfO( fr.) ^ pressure. For simplicity, we take A( f) 
as unity, and a pressure sensor can be simulated by inserting a suitable filter into 
the block diagram of Fig. 4. So, the output electrical signal will be numerically 
equal to the velocity potential. The receiver array, like the transmitter, is a 
rectangular array of point sensors of dimension M, x N, , parallel to the XY plane, 


with interelement spacings dyr and dy, meters, centered at rr = (2, Yr, zr). The 
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output electrical signal at each element (m,n), given by Eq. (37) with A(f)=1 


and r, = Enn = (mM dyr + Tr, n dyr + Yr, zr), iS 
Yan (J) = (finn) > (38) 
ery) =o (ee (39) 


d. Overall Transfer Function 
The overall transfer function relates the output electrical signal y, (t) 
at each point sensor (m,n) of the receiver array to the electrical input signal z(t) at 
the transmitter array. Substituting Eq. (39) into Eq. (32) we obtain, for the output 


signal, 


Yan (! =) Jo (LY) Hy (LYK aq) e] AU iv mn) q, df, (40a) 
uq = [xf ota Ey (SU Tm) p7] ZTV Inn E J?nf tag. (40b) 


The last expression can be rewritten as 


"OR [xo H( fr.) e 277 "af, (40c) 


00 


Or 


Yoon (t) = Fy pt XCF) AC Stan) } (40d) 


where H(fr.n) is the overall transfer function, given by 


H( fran) -f DE E Ean 170 Tong. (41) 


00 


Note that Eq. (40c) is valid only when the input electrical signal is independent of 
position and the receiver is a point sensor with unit frequency response, that is, 
A(f)= 1. On the other hand, Eq. (32) assumes only the independence of the input 


electrical signal with respect to position. 


B. THE RANGE INDEPENDENT ACOUSTIC CHANNEL 

In this section we will derive the expression for the overall transfer function for 
the particular case of a range independent medium. With reference to Fig. 3, the 
range independent medium characteristics change only with the Y (depth) 
direction, and the boundary conditions are independent of r and z, that 1s, the 
problem presents cylindrical symmetry. 

1. Medium Transfer Function 

In order to obtain the medium transfer function, we will solve the wave 


equation using a time-harmonic point source, that is, 


! 1 | 
V^p(4r) - Es p(tr) = ATE 6(r-15) , (12) 
c(y) dt 


zl 


with suitable boundary conditions . The solution, as shown in Eq. (6d), is of the 


form o(tr) = cies J Tht with &(r) satisfying the Helmholtz wave equation 
2 2 _ 
V^£(r) - K (y) €(r) 2 é(r-) , (43) 


where k(y) = 2af/c(y) is the depth dependent wave number. Applying the spatial 


Fourier transform with respect to r and z to Eq. (43), we obtain 
) Eg f + P zo) = A n + ha) Syy) (45) 
ay 


n SOIA 2) ert "ER oo €— Né ; , 
where kKy(y) 2 k(y)-4r (Kk + £^) and Z(f uw) 1s the spatial Fourier transform 
of E(r) with respect to 7 and z. Solutions to Eq. (44) are usually obtained by 
approximate methods, most notably the WKB approximation. We can write the 


solution to Eq. (44), without loss of generality, as 


(ky h) = dy) (A 39; (9) + B JO av + fis) (45) 
EX os) 
y 


The (transformed) velocity potential, obtained from Eq. (45) and the assumed 


(tr) 066) i rbt is given by 


(6f. y. E) = toy ty dm. (46) 


99 


The velocity potential y(tr) is given by the inverse Fourier transform of Eq. (46) 


with respect to rand z, namely 
/ n Pon f (fot) PAI, de. (47) 


Next, we will relate these results to those obtained in Section II-A -1c, the 
response of a linear time-invariant, space-variant system to a time-harmonic point 


source, as given by Eq. (6d), which can be written as 


ex^ 


T 


sn again PI il, qu 


00 


where fo = (i, yo, 2) = (5-25, y-Ys, 2-25) and dv=dk df, df. In order to 
compare Eqs. (47) and (45), we rewrite Eq. (47) using the term W defined in 


Eq. (45) as follows: 


E me ps 2 TL Is + L 25) ¿J27L 1 I Thak df, . (49a) 


"T _ li | y J 2T h Te RU E E de i (49b) 


Comparing the Fourier transform expressions in Eqs. (48) and (49b), the term Y. 


obtained from the solution of the inhomogeneous wave equation when the source is a 


unit amplitude, time-harmonic point source, is seen to be related to the medium 


transfer function through the transform 


V(ffoso uy) 7 8, el HG y, (50a) 


OI 
V (f f yo, fiy) - Fy Sho fy B39) e 3 27h Iodf, (50b) 


where y has been substituted for r as an argument of Ei in order to stress the 
depth-only dependence of the medium transfer function. Notice also that the 
arguments of V have been written in accordance with the linear systems notation 
introduced in Section II-A-1. 

The transfer function can then be obtained as the direct spatial Fourier 
transform of V( f fK.y;.f:y) with respect to. y, — y-y,. In the next section we will 
see that this last step 1s not necessary, so long as the transmitter array 1s composed 
of point sources. 

2. Output Electrical Signal 

The overall transfer function for a range independent medium can be 
computed from Eq.(41) by substitution of the expressions for D. (fv) and the 
medium transfer function. As we just mentioned, it will not be necessary to 
compute the medium transfer function if the transmitter elements are point sources. 
The directivity function for the transmitter array is a sum of terms of the form 


‘9 
2rv-I d 
Cat) e Pq, where Eng = (z, , Ya» is) ds the position of each transmitter 


element. Substituting that general summation term into Eq. (41) and expanding 
the exponential factors yields the result that the overall transfer function is a 


summation—over p and g—of terms of the form 


(D. Pre eo H (Shoh kiyn) e I Wn ugs ez dy 


OT 


MES cac E: 2 E s 
pal] CJ Tp) A E E A H hho hhi) d. 


“oa 


where r,2 n«md, and y,- y-4n dr are the (z,y) coordinates of the receiver 
array elements, and the (z,y) coordinates of the transmitter array elements are 
T= I +p d and y, = w*qd, . Recall from Section II-A-1 that the parameter 
To is the vectorial difference between observation point and source position. The 
integral in this last expression is the inverse spatial Fourier transform of the transfer 
function Fh fos A: Yn) with respect to £, f and f which yields a function of 
Ez). where L= Im=T, . Vo Un % > and 2 = z-z. But he transform 
with respect to f, is already given by Eq. (50b) as the function V(f fy f: y) 


obtained from the solution of the wave equation. Asa result, we can rewrite the 


above last expression as 


"pa | f Ylh yokiy) eI To eJ Thodi df. 


to 
cyt 


Therefore, the overall transfer function for the range-independent medium is given 


by 


Y; -1 Ni -1 
=> E | 
Li RE > > cool) i V (f fo go fit) € Th 
Ml Ni-1 OS 


pS eee | 
«cJ? lA. df d£. (51) 


The output electrical signal at each receiver position is given by Eq. (40c), 


repeated here for convenience: 


"NC =f XD Hr) °F fap, (40c) 


ec 


where H(fr,,,) is given by Eq. (51). Again, we stress that Eq. (40c) incorporates 


three assumptions regarding the linear system: 


e The input electrical signal is independent of position r, that is, z(&r) 2 z(t). 


e The receiver array elements are point sensors, that 1s, A, (51) = AE 


e The sensors have unit frequency response, that is, A(f)=1. 
Equation (51) assumes, in addition: 


e The transmitter array elements are point sources. 
e The medium has range independent characteristics and boundary conditions. 


e Transmitter array geometry is as described in Section I-A-2b. 
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a. An Alternative Representation 

We will now compare our results, as given by Eqs. (51) and (40c), with 
the results obtained by solving the wave equation for a range-independent medium 
in cylindrical coordinates and written in terms of a Fourier-Bessel integral. 
DiNapoli and Deavenport [Ref. 4:p. 80], for example, use that representation as a 
starting point for the description of numerical methods applied to range independent 
problems. 

The integral in Eq. (51) can be seen as the overall transfer function due 
to a single point source. From this point of view, H(fr,,) is just the weighted sum 
of the transfer functions due to each transmitter array element, as one would expect 
for a parallel connection of linear systems. These elementary overall transfer 


functions Ho (fr mn) are given by 
o oc E o 25) : 
“pa (ban - f f (f ho Yoki Yn) € pe € To q, d f. i (95) 


The overall transfer function can then be written as 


Mi -1 Ni -1 
= 


a 


Hftm)= y, A (53) 
ES Mi -1 n Nt -1 
D CS —- 5 


Let us make a change of variables in the integral in Eq. (52), exploiting 

the fact that the dependence of V(f£,vw,£;yn) with £ and f£ occurs due to the 
í ' í f 4) ") 

factor Ky) B Elo AE Rs in Eq. (44), that is, V is a function of (K^ - £^). 


to 
=~] 


Introduce the variables f. and y such that 
k = cosy and {= SUE (54) 


which is a change from "rectangular" coordinates (£,f) to the "polar" coordinates 


(4,7). With this change of variables, Eq. (52) becomes 


oc 


2T 
V ( f f yoiyn) | 25d £p roe c oS NEN Vay df, (55) 
0 


AO | 


0 


where 


V(ff.uow) 9 VG fu kiw! | (56) 
= o 
pod 
Equation (55) can be written as 
E Ea | 
im) = f W(ff.Yo:4n) | "e CIA suyo) y df ; (57) 
0 0 


where sina h ada hn and Rc 7 + 2°. Because of the 
periodicity of the sine function, the term a can be dropped from the exponent 
without altering the result. The last integral can then be recognized as the 


zero-order Bessel function of the first kind [Ref. 6:p. 184], that is, 


2r 
TOS | E, (58) 
0 
which leads us to the Fourier -Bessel integral representation 


TT) = 3 V(ff.,yoyn) J (27 R,) É df - (59) 


0 


Equation (59) is the same as Eq. (3.1) in [Ref. 4] if we perform the 
trivial change of variables k, = 2af. There, V is interpreted as the depth dependent 
Green's function. Equation (59) is the starting point to several approaches to the 
range independent problem, both numerical and analytical [Ref. 4:pp. 80-134]. 
Using the Fourier-Bessel integral representatiuü, the overall transfer function 


H( fran! can be written, from Egs. (53) and (59). as 


Mr) > y 2 zi Eg Yo Y) J,(OxER) f df. (60) 
| 0 


C. THE LAYERED WAVEGUIDE 
In this section we will derive the function Y for the case of a three fluid medium 


waveguide, commonly referred to as the Pekeris waveguide. 


1. Statement of The Problem 

Figure 5 shows the relevant physical aspects for our derivation. The point 
source is at position (z,,y,,z,) in medium II and the receiver is at (z,y,z) in the same 
medium. Medium II, with characteristic impedance pyc, is bounded at y= 0, the 
surface, by a fluid with characteristic impedance p,c, and at y= D, the bottom, by 
another fluid with characteristic impedance pc}. As described in Section B-1, in 
order to obtain V, we need to solve the Helmholtz wave equation, Eq. (44), and 
isolate the terms with y dependence from the solution, as shown in Eq. (45). An 
inspection of those two expressions show us that by setting z, and z, to zero, the 
solution of the wave equation will be automatically the function Y. The wave 


equations that describe the propagation in the three media are: 


O) 
2 v e Sv, e 0, (61) 
dy” 
9) 
Kyo Va + E Vo = OC y-ys) , (62) 
y 


and 


O 
e) ant 
ky Vy + yy =0, (63) 
dy 
2 DEO Ur ND ; 
where Ey; 2 kj-47 (bh + 6), k =27f/c; and the subscripts 1, 2, and 3 refer to 
the three different media. The argument of V has been dropped for simplicity. We 
seek the solution V, because we shall be concerned only with those problems where 


both source and receiver are in medium II. The boundary conditions are the 
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Fig. 5. The Three Layer Waveguide. 
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continuity of the acoustic pressure and the normal component of the acoustic 
particle velocity at each boundary. Acoustic pressure is related to the velocity 


potential by 


Pi(f) 5 32v fp; Vi, (64) 


and the normal component of the acoustic particle velocity at the boundaries is 


related to velocity potential by 


ely (65) 
dy 


2. Solution 
With the source located in medium II, the solutions in media I and III must 
be in the form of propagating waves moving outward from the boundaries. In 
media II, there are waves propagating in both directions—increasing and decreasing 


y—due to the reflections at the boundaries. The expected form of the solutions are 


Y= A, Ody <0, (66) 
Voa = Aga ez : B Boa e Ha E 0 $ y < Ys > (67) 
V op 7 Asp My "+ Bop e Ra € UNUM (68) 


and 


Vea ale (69) 


The solution in medium II has been split in two in order to account for the point 
source, as usual in the Green's Function determination. Two more boundary 


conditions are necessary: 


e Continuity of the solution V, at y= y,. 


e Discontinuity of the first derivative of V, at y = y,- 


The value of the discontinuity can be found by integrating Eq. (62) over a region 
| y-v,| « cand taking the limit as «> O. Performing this operation and taking into 


account that, in the limit, the integration of the continuous function V, is zero, we 








obtain 
iU cC d 
lim —s Vo dy=1, (70a) 
e+0 dy” 
Ug S 
OI 
d d 
D 2b UE d 2a Ye ( ) 
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a. Boundary Condition at y — 0 
In order to obtain the expression that represents the continuity of 
pressure at y = 0, we equate the value of pressures due to V, and V,,. Substituting 


the assumed solutions—kEas. (66) and (67)—into Eq. (64) yields 


Py A, m (Asa + Boa) i (71) 


For the normal component of particle velocity, Eq. (65) and the assumed solutions 


yield 


? A 


? ! wi 


Substituting Eq. (71) into (72). A, is eliminated and the relationship between the 


amplitudes of the solution Ya, 1s found to be 





Pi 
» pato o Mm 
7 E (n 
2a Po ko zr kyi 


When 4,» is real, corresponding to propagating waves in medium II, the factor Ro, is 
physically interpreted as the plane wave reflection coefficient at the surface. 
b. Boundary Condition at y = D 
Following the steps of the previous section, using the forms of the 


assumed solutions for V, and V,, at y = D, we find 
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pa Às RIT Po (Asp geo + By, e, (74) 
and 
-k,3 Az ete” = kyo (Aob RTN ¿ID (75) 


Substituting Eq. (74) into (75), 4, is eliminated and the relationship between the 


other two amplitude factors is found to be 


P3 





cmm 
Ao, OL Po y2 y3 = 
= Rag € IU zd 2 UR (76) 
Bop E k ES k 
p, Y? — v3 


When Kj is real, Ro, has the physical meaning of the plane wave reflection 
coefficient at the bottom. 
c. Source Condition 
Using the results of the two last sections, Eqs. (73) and (76). the 
solutions in medium II become 


oo (Myo + Ray Ayo! JU ELS Urs (0) 


a 


and 


Yop = Bob (Ass ¿ID JoY ats e Îy29 D Us < y < dela: (78) 
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The condition of continuity of V, at y — y, , applied to the above two equations, 


gives the relation between the amplitudes Ay, and By. 


Aga Fog ei 2ky3D 9hy2 Is 4 ea Ys 
= CA (79) 
Boy kyo Y 4 Ree -Jk yo Ys 





Finally, apply the condition on the derivative of V, at y- y, , as given by 
Eq. (70b), to obtain: 


Bae i E E ikg D E y2% Mya )- As. (e nsn, Ayo j == T 


which, together with Eq. (79), forms the necessary system of equations to obtain the 
two amplitude factors. B,, and 4,,. Upon substitution of the expressions for Bap 
and A,, into Egs. (77) and (73), we obtain the solutions V4, and V4, , which form 
the transformed transfer function — V(ff,w,f;v) we seek. Written in an 


abbreviated form, the final expression is 





/ 1 Pa 1 2k,oD QUE, y» -+ ¿Ho y> 
YU bode: fy) = 177 BOLD 
2 y? 1 _ og Ros € J + y2 
x (e 2, Re Jh yy (81) 


where OS y « D, O& y, € D , and, as usual, y> is the greater and y< the smaller 


between y and y, . 





3. Special Cases 
In order to check the solution for the layered waveguide, we will now 
derive the expression for the output electrical signal in two special cases whose 
analytical solutions are known. These results will be used to check the computer 
implementation. In both cases, the input will be a time-harmonic point source. 
a. Unbounded Homogeneous Medium ( Free Space) 
An unbounded homogeneous medium can be represented by Fig. 5 if 
the sound speeds c; and densities p; are equal for the three media. In this case, the 
reflection coefficients, as given by Eqs. (73) and (76). become zero. The 


transformed transfer function—EdQq. (81)——reduces to 





V( fb Yohiy) = i ey» ye (82a) 

ek yo 
ELA vob) = jd Apa ~ 9S) | (82b) 

ge 

ye 

OT 
re. o ES n 
Vf fgg) = j Prl al = ¡Lal (52) 
di ay 


Note that the right hand side of Eq. (82c) is expressed in terms of y, only, meaning 
that the unbounded homogeneous medium is space-invariant, as expected. 

The overall transfer function H(fr) for the case of a point source 1s 
given by Eq. (51) with M = M = c,,(f)=1. Substituting Eq. (82c) into Eq. (51) 


yields 


Ho = l j { jl old GU. id dj. (83) 
Figo OO 2k,9 
The above integral is the expansion of a spherical wave into plane waves [Ref. 7:p. 


B -5] and the resultant overall transfer function is 


-jk, R 
ATR 


H(fr) » -* 





| (84) 


where R= | r | = y i uA + E 


The output electrical signal y(t)—or velocity potential p(tr), in our 





(z-2,) yy) Haz) 
case—is given by Eq. (40c). Substituting Eq. (84) and X(f) = 6(f-6)—for the 


time-harmonic source—into Eq. (49c) yields 


D. ow | 
y(t) = e IT! = kB) 


(55) 
4r 
where k, 2 2zf, /c,. Equation (85) 1s just the expression for a spherical wave, since 
Has tlie distance between source and receiver. 
b. Surface Reflection 
Consider the case where media II and III are the same, that is, c = c 


and po = py. Asa result, Ro, is zero and V reduces to 


Vf fy 90,659) = ioe eRya9> (hy29< y po IZ) (86a) 
9 


€ y2 


Ce) imm (IYD, py al € Y) ). — (geb) 
y2 


OT 


Vlhko tokig) = j (ealt l a By Flyt 81) 07, 0%. (860) 
7 


Comparing Eqs. (82c) and (86c), we see that this problem can be treated by 
considering a homogeneous medium with two point sources, one located at y, with 
unit strength, and the second located at -y with strength Ra. The output 
electrical signal due to the first term in Eq. (86c) has already been computed in the 
previous section, and is given by Eq. (85). In the second term, we have the factor 
Ro, also a function of f, and f, through the dependence on k, and ko When 
performing the inverse spatial Fourier transform to obtain H( fr), it would be easier 
if R5, were constant. If we consider the boundary I-II as pressure release , that is. 


pulp = 0, then, by Eq. (13), Roy = -1. In this case, Eq. (86c) reduces to 


-Jka uc e. 
Y (ff Yo kiy) = j— (e J sale J y2l Yol) , (dI) 


where yi = y-(-y,) . From the results in the previous section, the output electrical 


signal is, by superposition. 


1 ; l . Ri 
(Ò = -— ROTA = E (88) 
K aA 
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where R' = || r-r! || = yz + yl 9 z; 9 Vn) (y 9) PUN: . If the distance 
R is large compared to the depth of the source y., a simplification is possible. From 
Fig. 6, for large R, Rz R,-AR, R'sR, « AR and ARz y, sinü— y y / BR, 
where R is the distance from the receiver to the midpoint between the source and 
in the amplitude 


its image. With these approximations and letting Rz R'a R 
factor, Eq. (88) becomes [Ref. 7:pp. 170, 408] 


S 


1 ko Ye Yr 
y(t)= -J RUOLI - Kf) sinl —— |. (89) 
PE dii R 


S S 


D. SUMMARY 
The transformed transfer function V(f£f,Yo,f£;y) for a range-independent 


medium is the solution to the inhomogeneous Helmholtz wave equation 


9 de 
ky(y) Y + d Y = d(y-ys), (90) 
y 


with suitable boundary conditions, where E = 4 ais £^ and 62x49 
Under the conditions stated in Section II-B-2, related to Eq. (51), the overall 


transfer function for a range-independent medium is given by 


M,-1 N.-1 
00 00 
-79 
Hl tmn) = x è ss f Í V(f foo sg) ITA 
Met NA Voto 
pS ee 


x eI Thode df (51) 
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Fig. 6. Geometry of Surface Reflection. 
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where the transmitter array is rectangular, centered at (z,,y.,z,), parallel to the XY 
plane, with odd dimensions Mt x Ni , complex weights wae and with 


elements positioned at Tp = Nn*tpr, Yq = Yet WY » ane 92 = ca 


Do Also. Tr inn 


m py 
Yo = YW-Yy, and Z% = 2,-2,. The overall transfer function H(fr,,,) can also be 


expressed as 





M-i N -1 
00 
Hfr)=2 Y Y a0 [ vasa) 
EMI NS: 
e n 
x J. (2r R.) É df, (60) 
where f 2. f^ + [is y E = 
VSL Yon) = Y(Sh Yo fin) : (56) 
AU a 
m 





The output electrical signal y(t), under the same conditions stated for the 


overall transfer function, 1s given by 


WE Í XD) Hle) I2 ay, (406) 


oc 


where H(fr,,,) is given by Eqs. (51) or (60). 


The layered waveguide illustrated in Fig. 5 is characterized by the transformed 


medium transfer function 


21 fa E guy Maio 
Wood) e je > 
y? IE Hone ila 


(Pp, Rca), (81) 


which is valid for 0 « y € D, 0€ y, € D, and where R,, is given by Ea. (73) and Ros 
by Eq. (76). If all three media have the same acoustic characteristics, then Eq. (81) 


reduces to 


. 1 E fc ‘ 
iGo = rl yal Yo | (82c) 


“= ve 
and the corresponding output electrical signal 1s given by 


1 (O = 
A - me AS Ro) (85) 
iT 


where R is the distance from source to receiver. In the case where media II and III 


are the same and the surface is pressure release, the function V is given by 


Te -Jk ol Y 
Veto) m $5 — (Cl ate. (87) 
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where y} = y+ y,- The output electrical signal, when the source -receiver distance 


is large compared to the source depth, is given by 


l y ? ko Ys Yr 
y(t) = -j— ROLIN ko Ro) sin| ———|, (89) 
2xR R 


S S 


where R, 1s the distance from the receiver to the midpoint between the source and 


its surface image. 
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HI. IMPLEMENTATION AND RESULTS 


In this chapter, a computer implementation of the range independent medium 
equations is described. Although only the medium transfer function of the layered 
waveguide is implemented, the modular nature of the linear systems approach 
enables us to analyze the results and draw conclusions valid for the more general 
case. Equation (60) illustrates this point. The summations and the coefficients 
pem that equation represent the transmitter array. The Bessel function 
J (2xf R,) is the range dependent factor. The function V(f f, y. f; y)—írom now 
on referred to as the medium transfer function, for simplicity—is the (range 
independent) medium dependent factor. 

Our purpose is to have a working algorithm to test the results from the previous 
chapter. Care has been taken to enable the implementation of other medium 
transfer functions, allowing for depth dependent speed of sound, so that the program 
can be used as an academic tool. Speed and style were secondary objectives. 


The implementation is based on Eq. (60), the Fourier-Bessel integral 


representation. 


A. IMPLEMENTED EQUATIONS 
1. Medium Transfer Function 


The medium transfer function ¥(ff,y,;y) is an oscillatory or exponential 





function of k,,—as seen in Eq. (81) —which, in turn, is a function of TN £^ - 


D : : 
= f. The function V is easier analyzed in terms of Kyo than in terms of f = 
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= (ffe)*- (k,9/27)" , which is nonlinear. Therefore, V will be treated as a 


function of ko, and the integral term in the expression for the overall transfer 


y? 
function, Eq. (60), will be modified accordingly. 
a. WV as a Function ofk, 
The function V is already expressed in terms of kyo in Eq. (81). The 


reflection coefficients, on the other hand, are expressed in terms of kr k 


yn and k,3: 


In this section we establish the relationships between Kyo and the quantities f, k 
and ka . These relationships will enable us to represent the reflection coefficients in 


terms of the variable k,;, as well as to change the variable of integration in the 


y? 
expression for the overall transfer function to hyo: 


The correspondence bet ween Kyo and f is one-to-one. When f < (f/a) , 


kyo is taken as the positive root 27 (f/c) E f^ , a consequence of our choice of 
solutions, Eqs. (67) and (68). When f > (f/c,), k,o is imaginary, corresponding to 
an exponentially decaying V. When all factors in Eq. (81) are multiplied out, V is 


seen to be composed of a sum of factors of the form 


-K aa 
A e ya 


where a is a nonnegative real number. In order to have V be a decreasing 
exponential, the negative root must be chosen when ko Is imaginary. Similar 
reasoning yields the same results regarding k,, and k,3) which can be summarized as 


follows: 
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Lu. 


y (91) 


+ EA ok Sek. 
-;| 2-2 OR ki 
where t= 1, 2013, hk, =2xf and k; = 2zf/c,. Using Eq. (91), the wavenumber 


y-components ks and k,3 can be expressed in terms of kyo as 


| Dan? 2 

| o y Ko K<k 

kyp = | (92) 
D 


: | 12 2 2 
= Boa a r p 


where p=1lor3, and k = | ks — E . Using Eq. (91), the integrations over f 
can be transformed into integrations over ko. Using Eq.(92), the reflection 
coefficients R,, and R,, . as given by Eqs. (73) and (76), can be expressed in terms 
of ko. 
b. Air, Water, Fast Bottom Waveguide 

The present implementation of the medium transfer function assumes 
that cj< co. When the surface is almost pressure release, which is the case of the 
air-water interface, and the bottom is fast, c> c , the medium transfer function is 
composed of many sharp peaks, which presents some difficulties for integration 
routines. 

(1) Poles of the Medium Transfer Function. When the surface is 


pressure release, that is, Ay, = -1 , Eq. (81) reduces to 


ER ed 2kya D NL Aya | 
V(fkoyoy) T - RS AAA sin( ko y<), (93) 
kyo i+ Ras e A 


which is valid for O « y 4 D and O& y, & D. If the bottom is fast, k > kz, and, 


from Eq. (91), there is a range of values of k, for which k,, is real while k,3 is 


imaginary, that is, kg < kp < kọ or O< kyo <y hy - ky. For this range of values 


of kyo, the reflection coefficient Rog, given by 


P3 
p A Ore ks 
Ros = que > (94) 
Da kyo + kg 
becomes a unit magnitude complex number, 
Ry, = 0%, (95) 
where 
2 2 E 
pa Im( ky j Pa | kg — Ky — hs 
ino dee > | (96) 
P3 Kyo P3 kyo 


In this case, the denominator of the expression for V, Eq. (93), may have zeros, 


given by 


D 


1 + Bos cJ? yo umm (97a) 
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1467 XD - 9) «9. (97b) 





Therefore, 
¢ - kD = (2n+ 1) e n=0,1, 2... (98a) 
or 
tan kD = » n n — (98b) 
| EEE 


3 
T 
2 
> 
$ 
j 
4 
2 
3 
3 
> 
L 
3 


P3 0 
tan à zm -— ————, (99) 


P» DE 


where 0= koD and a= " K - E . Solutions to the above transcendental 
equation in terms of k,, = 6/D are the poles of the medium transfer function, for 
the case of a pressure release surface, Eq. (93), and a fast bottom. The physical 
interpretation of these poles is that they represent the discrete eigenvalues kan 
corresponding to the trapped or normal modes in the waveguide[Ref. 8:pp. 430-440]. 
The graphical solution of Eq. (99) [Ref. &:p. 438] reveals that: 
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e Poles exist only if a> 7/2. If a} x , the poles will be separated by Ak. o= 7/D. 


e The total number of poles is given by the integer number equal to [(a/7)- 0.5] 
or, if the result is not integer, its nearest larger integer. 


When the surface is almost pressure release, that is, Ry, 7-1 , 
the medium transfer function presents peaks at values of k corresponding to the 
poles computed as above. 

2. Overall Transfer Function 
In order to change the variable of integration in Eq. (60), the interval of 
integration f. € (0,00) is subdivided into (0, f/c)) and ( f/c,, oo ). After the 


change of variables according to Eq. (91) and a rearrangement, Eq. (60) becomes 
zc EN TUA CUP > (100 
Dinn Da = T 4 i V hA ym o Sn Vo FETO kyo dk» ; (100) 


where the path of integration I’ consists of the imaginary axis over the interval 
(¡o , 30) plus the interval (0, f/c) on the real axis. The summations over p and 


g are as indicated in Eq. (60). 


B. THE SUBROUTINES 
In this section we describe the subroutines and present the results for some 
examples. Also, a description of the main program and some auxiliary subroutines 


is given, which will help to clarify some aspects of the implementation!. 


l'l'he main program and auxiliary routines were written by Dr. Lawrence J. 
Ziomek, at the Naval Postgraduate School, as part of an ongoing research project in 
model-based signal processing. 


1. Main Program and Auxiliary Routines 

The main program's simplified pseudo-code is shown in Fig. 7. It consists 
basically of calls to subroutines. The readin subroutine reads the input data, 
consisting of transmitted signal characteristics, transmitter and receiver array data, 
medium characteristics, and control data, including the ocean medium transfer 
function (OCNTF) variable used to select a transfer function, when more than one 
type of transfer function is implemented. The ceq subroutine generates the input 
electrical signal z(t). The chfmn subroutine computes the overall transfer function, 
Eq. (100), and its estimated error. The phjmn subroutine provides a tabular output 
for the values of the overall transfer function and its error. The calyce subroutine 


computes the output electrical signal, according to Eq. (40c). 


Main Program  /computes output waveform of an acoustic channel/ 
EXECUTE Readin  /read input data/ 

EXECUTE Ccq  /generate input electrical signal/ 

EXECUTE Chfmn  /compute overall transfer function/ 

IF (flag print is TRUE) THEN 

EXECUTE Phfmn /print overall transfer function/ 
END IF 
EXECUTE Calyce /compute output electrical signal/ 


End Main Program 


Fig. 7. Simplified Pseudo-Code of the Main Program. 
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a. Subroutine Ccq 

The subroutine ccg generates the input electrical signal z(t), using the 
complex envelope representation. From the specifications read as input data, this 
subroutine computes the time samples of a reference complex envelope and its 
Fourier series coefficients. The Fourier series is truncated, yielding the complex 
envelope of the transmitted signal. 

(1) Complex Signal Representation. An amplitude and angle 
modulated carrier r(t) can be represented by its complex envelope T(t) as [Ref. 


Sep AR 
at Re E (101) 
where 
0) = a(o) AB, (102) 


a(t) is the (real) amplitude modulating signal (Amp), 9(t) is the angle modulating 
signal (rad), and f is the carrier frequency (Hz). The Fourier transform of the 


signal is related to its complex envelope transform by [Ref 5:p. 187] 


R(f) = O.5{R(f- £) - RF'C-f- (103) 
For a rectangular -envelope LFM pulse, a(t) and 9(t) are given by 


A ,|t] < TP/2 


all) = ; 104 
M) ie AN dc 


02 


and 


T SW 
0(t) 2 ——— (2 = Dp A (105) 
TP 


where A is the (constant) amplitude (Amp), TP is the pulse length (s), Sw is the 
swept bandwidth (Hz) and Dp is the phase-deviation constant (rad s Both Sw 
and Dp are allowed to be either positive or negative, corresponding to the up chirp 
and down chirp LFM pulses, respectively. If the modulated signal r(t) is made 
periodic with period Ts, then the complex envelope can be represented by the 


Fourier series 


, 
AGE > by SAI (106) 
k= -f 
where 
1 qto! | 
b = -=f O N (107) 
oo hee 


and f,=1/T,. Equation (106) assumes a complex envelope bandlimited to Af. 


From Equation (106), the transform of the complex envelope is given by 


RU) V, hk. (108) 
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If the complex envelope is sampled with sampling period Ts = To/(2K + 1), which 
is basically the Nyquist rate, then, from Eqs. (106) and (107), 


F(T.) = Y RE (109) 
= 


and 


f 
1 n 
E È | A pr 2xkl/ L . (110) 


where 7, = T./L is the sampling period, and £=2K+1 is the total (odd) 
number of time and frequency samples. From Eqs. (103) and (108), the signal 


transform is given by 
k 
R()=05 5 (dbf REY + BE SCSI, (111) 
ieri 
from which the sampled signal can be derived as 
A a a 
ye Rel Y 4, J TkL Jrj | (112) 


TE, 


(2) Subroutine Description. Figure 8 lists the pseudo-code of the 


subroutine ccg. The main input data for the subroutine are: 


Pulse length TP, and pulse period (or pulse repetition period) To. 
Swept bandwidth Sw and pulse amplitude A. 


Maximum order of frequency component to be processed, Kmax. The number of 
frequency components to be processed is (2 Kmax + 1). 


Sampling rate multiplier SR, a term to be multiplied by the Nyquist rate in 
order to obtain the actual sampling rate. 


Carrier frequency f. . 


A flag test used to indicate option for plotting. 


Subroutine Ccq 

Compute data for reference signal: 
total number of samples L for the complex envelope /Eq. (114)/ 
sample period TS for the complex envelope /Eq. (113)/ 
phase—deviation constant DP /Eq. (105)/ 

Generate samples of ref. signal complex envelope /Eq. (102)/ 

Compute complex Fourier coefficients /Eq. (110)/ 

IF (flag test is TRUE) THEN 
Compute and plot transmitted signal /Eq. (116)/ 

END IF 


End Subroutine Ccq 


Fig. 8. Simplified Pseudo-Code for the Signal Generator Subroutine. 
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The total number of samples for the complex envelope is given by 


L= To/Ts, (113) 


where To and T, are the pulse repetition period and complex envelope sampling 
period, respectively. The sampling frequency f =1/T; is the product of the 
sampling rate multiplier SR by the Nyquist rate, twice the complex envelope 
bandwidth. In terms of the input parameters listed earlier, Eq. (113) can be 


rewritten as 


Ís 


ť— MÀ 
La Tort 29998 94] ZE (114) 


x bk(envelope) 


which is the expression used in the subroutine. If the number of samples is less than 
the number of frequency components to be processed, (2 Kmax + 1), £ is increased 
accordingly. If the computed value of L is even, the result is increased by one. The 
phase -deviation constant DP and the sampling period Ts are computed according to 
Eqs. (105) and (113), respectively. 

The subroutine proceeds to compute the samples of the complex 
envelope of the reference signal —using Eqs. (102), (104) and (105)—at the time 
instants t= | Ts. The Fourier coefficients are computed next, using Eq. (110). 

The signal to be transmitted is defined by multiplying the Fourier 


coefficient series {6,}—-computed for the reference signal—by a window of length 
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(2 Kmax + 1). The simple truncation, which corresponds to a rectangular window, is 
avoided because of the ringing it would cause on the waveform. Instead, the 


Hamming window 
w(n) = 0.54 + 0.46 cos(a nf Knax) , -KmaxS 2 < Kinax (T5) 


is used. Therefore, the transmitted signal and its complex envelope are given by 


mas 
E | y uh) bd STR E Jim -K«ICK, — (116) 
4 A max Á 


and 


Ímax 
HIT.) = Yuh A reter, aim 


ke -Ímax 
The transform of the transmitted signal 1s given, from Eq. (116), by 


Ámax 
XD =05 5 {MSC LA + BE SILA}. (118) 
k = -Å max 


While the sampling rate is high enough for the complex envelope 
representation, the signal x(t) usually requires a higher rate. In order to recover the 
signal from the Fourier coefficients, the series (5,) must be padded with zeros, so 


that the modulated signal samples occur more than twice per carrier period. Both 


DN 
msg 


the subroutines ccg and calyce, to be described next, use six samples per carrier 
period for plotting purposes, for a total of 6f T; samples on the interval 
(-To/2, To/2). This zero padding is accomplished by using the value Le = 6f T; 
instead of L in Eq. (116). 
b. Subroutine Calyce 

The subroutine calyce computes the output electrical signal, both in 
complex envelope y(t) and in band limited y(t) forms. 

(1) Sampled Output Signal. The electrical signal y,,(t) at the 


output of the element (m,n) in the receiver array is given by Eq. (40c). 


(40c) yields 


Es 
om emen] Y uk) b ir SL go 


Í max 


after performing the integration indicated in Eq. (40c), and where f T, — lfL. 
From the definition of the complex envelope, Eq. (101), the sampled complex 


envelope of the output signal is given, by inspection of Eq. (119), by 


x ] 
DS y wk) b, AC + kr) STE (120) 


k = -Imax 


Equations (119) and (120) are both implemented in the subroutine calyce The 
observations made about the recovery of the signal from its Fourier coefficients— 


regarding the zero padding of Fourier coefficients—are also valid for y(t). In order 
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to have a correct reconstruction of y(t), the constant L in Eq. (119) is substituted 
by Le = 6f. To, corresponding to six samples per carrier period. 

(2) Subroutine Description. Figure 9 lists the pseudo-code of the 
subroutine calyce. The main input data for this subroutine, besides those listed for 
the subroutine ccq, are a three-dimensional array with the values of the overall 
transfer function for each frequency (f. + kf) and position (m,n) in the receiver 
array, and the Fourier coefficients b, computed in the subroutine ccq. 

The computation of the total number of samples L and sampling 


period 7; are as described in the previous section, in connection to the generation of 


feio 


the transmitted signal. The complex envelope samples Y, (175) are compute 


according to Eq. (120). The signal samples y, (!T;) are computed using Eq. (119), 


with zero padding. 


Subroutine Calyce 
Compute: 
total number of samples L for the complex envelope 
sample period TS for the complex envelope 
Compute complex envelope samples  /Eq. (120)/ 
Compute and plot signal  /Eq. (119)/ 


End Subroutine Calyce 


Fig. 9. Simplified Pseudo-Code for the Subroutine Calyce. 
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2. Computing the Overall Transfer Function 

The computation of the overall transfer function is performed by a group of 
subprograms subordinated to the subroutine chfmn, called by the main program. 
The relationships between the subroutines used for this implementation are shown 
in Fig. 10. The program was implemented in FORTRAN, with most of the 
variables shared through common blocks. 

The pseudo-code for the subroutine chfmn is shown in Fig. 11. This 
subroutine implements Eq. (100), with the path of integration F truncated to the 
intervals (-jlower , 70) plus (0 , upper). 

3. Integrand Implementation 

With reference to Fig. 10, the subprograms related to the computation of 
the integrand are reinteg, iminteg, integr, refl, hwvgky and dbsj0!. They are all 
implemented as FORTRAN external functions. 

The functions reinteg and iminteg just return the real and imaginary parts 
of the integrand, respectively. They are called by the integration subroutine, whose 
first argument is the name of the function to be integrated. 

a. Function Integr 

The function integr computes the (complex) integrand in Eq. (100). 
The pseudo-code for this function is shown in Fig. 12. The argument E, is 
interpreted as imaginary when it has a negative value, in all subprograms. The IF 


statements marked in the list are related to the medium transfer function. New 


IIMSL, Inc.. DBSJO, SFUN/ Library, Bessel Functions of Integer Order, 1984. 
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Fig. 10. Structure of the Subroutine Chfmn. 
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Subroutine Chfmn 
FOR ( each frequency f= £f+X£) DO /-Kmax SKK Kmax/ 
EXECUTE Weight /compute transmitter array weights c__(f)/ 
FOR ( each position y, in receiver ) DO 
EXECUTE Break /compute breakpoints for integration/ 
Compute upper and lower limits of integration 
FOR ( each position x in receiver ) DO 
EXECUTE Qinteg(Real part) /compute Re{ H(f,rmn) } 
and its error Rel FHC f, rn) +/ 
EXECUTE Qinteg(Imag part) [compute Im( HO ERN 
and its error Im{ EHC f, Imn) }/ 
Compute magnitude and phase of Hand EH 
END FOR each zr, 
END FOR each y, 
Plot magnitude and phase of H versus position (r, y) 
END FOR each frequency 


End Subroutine Chfmn 


Fig. 11. Pseudo-Code for the Subroutine Chfmn. 


Function Integr(k,») 
/this function is executed for constant frequency, receiver 
position (m,n) and k,o/ 
Compute k, /Eq. (92)/ 
IF! (layered waveguide) THEN 


Ro, = Refl (surface data) and A.,= Refl (bottom data) 
END IF 
/compute product V Kk o. for each y,, as an array Ftran( q)/ 
FOR ( each position y in transmitter array ) DO 
IF! (layered waveguide) Ftran(g) = Hwvgky /kya WCS, koos Yoi Yn)” 
END FOR each fe 
/compute integrand according to Eq. (100)/ 
FOR (each position T, in transmitter array) DO 
Compute Bessel function J(k A.) /R, is the horizontal 
distance from transmitter element to receiver position/ 


FOR ( each position y, in transmitter array ) DO 


ácumulate product — c4(f)-PtranCq) - J, (E. R,) y 5 


al 
M 


END FOR each Ya 
END FOR each To 
Integr = (acumulated product )/27 


End Function Integr 


lOUther IF statements may be necessary for different medium transfer 
functions 


Fig. 12. Pseudo-Code for the Function /ntegr. 
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transfer functions are implemented by adding new IF statements. The functions 
integr and break, to be described later, are the only ones that need to be modified 
when incorporating other medium transfer functions to the program. 
b. Function Refl 
This function computes the complex reflection coefficients. Its 
pseudo-code is shown in Fig. 13. 
The expressions for both reflection coefficients, Eqs. (73) and (76), can 


be written as 


Function Refl(c,, > Posie kyo) 
/computes reflection coefficient A9p/ 
Compute k,n /Eq. (92)/ 
IF (the speeds of sound are different) THEN 
Compute reflection coefficient by Eq. (121) 
ELSE 
Compute reflection coefficient by Eq. (122) 


END IF 


End Function Refl 


Fig. 13. Pseudo-Code for the Function Refl. 
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where Kyp Is given by Eq. (92). If the speed of sound is the same in both media, 
then Eq. (121) can be simplified to 


Ry m. (122) 


c. Function Hwvgky 
The function hwvgky returns the complex product ko V( f ko Von)» 
the medium dependent factor in the integral expression for the overall transfer 
function, Eq.(100). This function is specific to the layered waveguide and 


implements Eq. (81), multiplied by k 


ya It is called by the function integr, as 


described earlier (see Figs. 10 and 12). Its implementation is trivial. Its arguments 
are the frequency, ko and the y coordinates of the transmitter and receiver array 
elements, y, and y, , respectively. The reflection coefficients Rọ, and Raz, and the 
depth D are available through common blocks. 
d. Examples 

In this section we examine the behavior of the integrand in the 
expression for the overall transfer function, Eq. (100). In all plots in which k,, is 
the independent variable. negative values of k, are imaginary and positive values 
are real. For simplicity, Sed negative imaginary and positive real values of ko are 
plotted on the same axis. The symbol "0'' appearing on the plots is used by the 
plotting routine to mark different plots on the same graph, but it has no special 
meaning in the present context. 

(1) The Range Factor J(k Ro). Figures 14 and 15 show the plots of 
J((KR,) vs k 


yo 


which is related to k through Eq. (91). The frequency is f= 
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Fig. 14. Range Factor J(k R.) for f= 500 Hz, q = 1,500 m/s, and A, = 100 m. 
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Fig. 15. Range Factor J(k R,) for f= 500 Hz, c, = 1500 m/s, and 
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= 500 Hz and the sound speed is c = 1,500 m/s, which results in a wavenumber 
k, = 2.09 m. In Fig. 14, the distance is R, = 100 m and, in Fig. 15, it 18 1,000 m. 
In those figures the range of kyo is (-70.25 &j , 0) plus (0,4). An important 
characteristic of the factor J (k R,) is that it oscillates slower near kyo = 0. 
Another characteristic behavior is that it oscillates faster as the range increases. If 
the other terms in the integrand also have "slow" oscillations, then the greater 
contribution for the result of the integration comes from the region kj, x0, as a 
result of the smoothing nature of the integration operation. Physically, the plane 
waves that travel near the horizontal, that 1s, those with Kyo x 0, contribute more 
to the total field than those that interact more with the boundaries. In order to 
compare the "frequency of oscillations" of the Bessel function and oí the medium 
transfer function, recall that, for kR,> 27, the zeros of J(kR,) occur 


approximately at intervals of 





bk = —. (123) 
Ro 
(2) The Medium Factor. Figure 16 shows the plot of the term k,, Y 
for the layered waveguide whose characteristics are shown in Table 1. The depth D 
of the ocean is 100 m and the frequency is 500 Hz. The position of the point source 
is (ts, Ys, 25) = (0, 30, 0) meters and the point receiver sensor is at (Tr ‚Yr. Zr) = 
= (1000 ,50 , 0) meters, corresponding to the range factor shown in Fig. 15. 

This waveguide is of the | almost -pressure release -surface/fast - 
bottom type. The peaks are evident in the figures. From the discussion in Section 
III-A-lb, there are 33 peaks in the range 0< Kyo < 1.04 m^! , separated by 
Ak, = 7/100. 
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Fig. 16. Waveguide Medium Transfer Function Times ko: (a) Real Part(cont.). 
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TABLE 1. WAVEGUIDE ACOUSTIC CHARACTERISTICS 


Medium Speed of Sound Density 
sii : 
(ms) (Kg m^) 
I 343.0 qe» 
I] 1500.0 1026.0 
III 1730.0 2070.0 


A criterion to set an upper limit of integration in Eq. (100) —the 
upper discussed in Section III -B-2——ould be based on the comparison between the 
"periodicities" of the medium and range factors in the integrand. If the factors in 
the expression for V, Eq. (81), are multiplied out, the faster oscillatory factor—the 
one with the greater exponent——has a "period" of z/ D or more, with respect to ka. 


that 1s, 





E 3 (124) 
D 
From Eq. (91), the increments ÓK e and ók, are approximately related by 
Kyo 


when k, is real. Equations (123) and (125) can be used to compute the interval 
between zeros of the range factor, when it is plotted versus ky». If we require, for 
example, that the integral be truncated when the Bessel function passes through 
zero N times during one "period" of the medium factor, then, substituting Eq. (123) 
and 0h, =7/ ND into Eq. (125), yields 


ky 
upper = ———————— . (126) 


/1 + (Ref NDS? 


The upper limit computed according to Eq. (126) can smaller than the last peak 


position, given by y ky - ks [see discussion following Eq. (99)]. In this case, the 
limit could be increased, or those peaks outside the limit could be simply discarded. 
In this implementation, the upper limit is set to ka. 

The lower limit of integration, in the region where Kk is 
imaginary negative and |X,9-Y | decreases exponentially as k, - y, should be 
set to minimize the error of integration for the slower decaying factor in ko- Y—the 
one with the exponent of smallest absolute value. The worst case is when there is 
one constant term, which happens when the source and the receiver are at the same 
depth. for example. Note that the term kyo is pure imaginary when Kyo 15 
negative imaginary, and decreases towards zero as LES ES e 20lkyal The lower 
limit is set to -70.125 k, in this implementation. 

(3) The Integrand. Figure 17 shows the plot of the integrand for the 
waveguide discussed above, at f= 500 Hz, for a range Ro = 1,000 m. For k 
greater than about 1.5, corresponding to N x 10 in Eq. (126), the integrand is fairly 


symmetrical around zero, suggesting a low value for the integral in that region. 
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17. Waveguide Integrand, Range Ro = 1,000 m. (a) Real Part. 


0.25 


0.20 


0.15 


REAL( INTEGRAND ) 


0.10 


0.05 





0.6 0.8 1 1227 1.4 1.6 1.8 2 gu 
KY2 (1/METER) 


Fig. 17. Waveguide Integrand, Range Ro = 1,000 m. (a) Real Part(cont.). 
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4. Other Subroutines 
a. Subroutine Break 

This subroutine divides the interval of integration in order to decrease 
the error of the integration subroutine. The criteria used in this division are the 
behavior of the medium transfer function and the relationship between the 
‘‘periodicities’’ of that function and of the Bessel function in the integrand. For 
each different medium transfer function, a different interval subdivision strategy 
must be implemented. The upper and lower limits of integration are independent of 
the type of medium, in this implementation. Nevertheless, they could be made 
medium dependent, in which case their computation should be included in the 
subroutine break. Figure 18 lists the pseudo-code for this subroutine. 

In the case of the lavered waveguide function, the breakpoints are set 
equal to the locations of the peaks. if any, as discussed in Section III-A -1b, 
regarding the almost - pressure-release-surface/fast -bottom. "The interval of values 
of k, outside the range of the peaks is divided evenly into subintervals of size 
ox / D, corresponding to five "periods," as given by Eq. (124). 

(1) Subroutine Description. The input parameters for the subroutine 
break are the frequency f the speed of sound c,—computed in the subroutine chfmn 
and taking into account the gradient of the speed of sound, for the case of other 
media. Other variables shared through common blocks are the type of medium, 
OCNTF, and the waveguide characteristics. The outputs are the number of 
breakpoints and a vector with the values of ky» at the breakpoints. 

Initially, the subroutine divides the interval of integration evenly. 
Then, it proceeds to check for the presence of peaks in the medium transfer 


function. as discussed in Section III-A-1b. The peaks are the zeros of Eq. (99), 
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Subroutine Break(/, c, , breakpoints , number. brkpts) 
IF! (layered waveguide) THEN 
[divide range evenly (-0.125 k, O)U(CO, k)/ 
number brkpts= (1.125 k, + A ko) -1 
FOR (n=1 to number brkpts) DO 
breakpoints(n) = (0.125 ko) + n Akys 
IF (breakpoints(n) < 0) inder = n [index keeps track of 
last nonpositive element of vector breakpoints/ 
END FOR n 
/check for presence of peaks/ 
IF (peaks do exist) THEN  /a» x22, see Eq. (99)/ 
/the peaks are located at k., 7 0: D, see Eq. (99)/ 
number peaks- [Ca[zy)- 0.5| / discussion after Eq. (99)/ 
FOR (n= inder+1 to index+1+ number_peaks) DO 
Compute n_th_zero of the function Denwgd 
breakpoints(n) = n_th_zero + D 
END FOR n 
Update number_brkpts to reflect total number of breakpoints 
END IF 
END IF 
End Subroutine Break 


lüther IF statements may be necessary for different medium transfer 
functions 


Fig. 18. Pseudo-Code for the Subroutine Break 


implemented as the FORTRAN external function denwgd. The zeros are found by 
the subroutine dzbren!. The (positive) breakpoints originally computed are changed 
to the values of the positions of the peaks. This implementation assumes a lower 
limit of integration of -70.125 k, and an upper limit equal to k,. 
b. Subroutine Weight 2 

The subroutine weight computes the transmitter array complex weights 
cpa(f). The input parameters are the frequency; the transmitter array data; the 
direction cosines u' (with respect to the X direction) and v' (with respect to the Y 
direction) of the transmitter array main lobe; the speed of sound at the center of the 
array—set to Cp for the layered waveguide, but included in the present 
implementation to allow for other media—and the gradient of the speed of 
sound——set to zero in this implementation. The complex weights are given by [Ref. 


Dap 13a} 


exse e (127) 

where 
apq = Ap Ba, (128) 
pa = E (ul pd v adj), (129) 


IIMSL, Inc., DZBREN, MATH] Library, Nonlinear Equations, 1985. 


“This implementation is based on a program originally written by Dr. 
Lawrence J. Ziomek, at the Naval Postgraduate School. 


Col Ya) 


pee 


| (130) 
f 


de and B, are the amplitude weights along the X and Y directions, respectively, À 
is the (depth-dependent) wavelength, d,, and dy, are the interelement spacing 
along the X and Y directions, respectively, and c,(y,) is the speed of sound at the 
depth Y, 

c. Subroutine Qinteg 

The subroutine ginteg is an interface for the subroutine dg2ags!, which 
actually performs the integration in Eq. (100). The input arguments are the name 
of the external FORTRAN function to be integrated, the limits of integration, the 
breakpoints, the number of breakpoints, and the required error in the result. The 
output arguments are the result of the integration and its estimated error. 

Initially, each interval of integration, defined by the breakpoints, is 
subdivided in order to improve the convergence of the subroutine dg2ags, which is 
called by ginteg to perform each intermediate integration. If the error of the 
integration of a subinterval, as computed by dg2ags, exceeds the requirements, that 
subinterval is further subdivided. Finally, the results and errors of integration 
provided by dg2ags are accumulated. 


As indicated in the pseudo-code of the subroutine chfmn, Fig. 11, each 


IMSL, Inc., DQ2AGS, MATH] Library, Integration and Differentiation, 
1985. 
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integration is performed twice, for the real and imaginary parts of the integrand. In 
terms of computing time, this is the most critical section of code. Not only does the 
subroutine ginteg calls dg2ags many times—at least ten times between breakpoints, 
in the present implementation—but, what is worse, dg2ags can execute the function 
to be integrated, ultimately integr, hundreds of times. The computation time 


increases with frequency and range: 


e The interval of integration increases as O(k,), directly proportional to 


frequency. 
e As the range increases, the range factor J,(k;R.) oscillates faster, degrading the 


convergence of the ite ahga subroutine, causing more executions of the 
function integr and more frequent interval subdivisions. 


The effect of range can be ininimized if the limits of integration, lower and upper, 


decrease (in absolute value) with range, as in Eq. (126) for the upper limit. 


C. RESULTS 

The algorithms described in the previous section are tested by using two 
examples for which the results are easily interpreted. The unbounded homogeneous 
medium was presented in Section II-C-3a as a particular case of the layered 
waveguide, all three media with the same acoustic characteristics. The surface 
reflection problem was analyzed in Section II-C-3b. In both cases, the input to the 
system is a time-harmonic point source and the output electrical signal— 
numerically equal to the velocity potential—is given by Eqs. (85) and (88). From 
Eq. (40c), the time-independent term of the output signal is, for a time-harmonic 


input, equal to the overall transfer function, that is, 
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"T [xo Hx, ,) e 271 a, (40c) 


00 


"E fura Hfr) It af, (131a) 


Yn C!) = HP Stahl, (131b) 


Therefore, the magnitude of the overall transfer function 1s equal to the magnitude 
of the velocity potential. 


Results for the waveguide example used so far (see Table 1) are also presented. 


Figure 19 presents the plot of the magnitude of the overall transfer function 
as a function of horizontal distance r from the source. The acoustic characteristics 
are those shown in Table 1 for medium II. The depth of the source is 30 m and of 
the receiver array, 50 m. The expected result is a straight line on a log-log plot, 
corresponding to a velocity potential proportional to the inverse of distance. At 100 
m, the error of the computed value is 6%, and increases with distance up to 29% at 
3,100 m. The integration routine computed the integrand 3,360 times (total for 
real and imaginary parts) for 100 m of distance, a number that increased to 34,776 
at 3,100 m, indicating an increased difficulty in convergence of the integral. An 
analogous computation for the layered waveguide took 45,000 and 70,000 


computations of the integrand, for 100 m and 3,100 m, respectively. 
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Fig. 19. Overall Transfer Function for a Homogeneous Medium. 


2. Surface Reflection 

Figure 20 shows the overall transfer function as a function of depth, for a 
distance of 1,000 m between source and receiver. The acoustic characteristics of the 
two media are those shown in Table 1 for media I and II. The surface is, for all 
practical purposes, pressure release. The depth of the source is 30 m. For these 
data, the approximations for Eq. (89) are valid, indicating a sinusoidal variation 
with depth. The error increases from zero at the surface, up to a maximum error at 
30 m, the depth of the source——see discussion in Section III-B-3d, regarding the 
lower limit of integration. At 30 m, near a theoretical maximum, the error is 4596 . 
At depths greater than 60 m, the error 1s negligible. For comparison purposes, the 
integrand was computed about 11,000 times for each value of the overall transfer 
function. 

3. Layered Waveguide: Waveform Prediction 

Figure 21 shows the transmitted signal. This signal was obtained by 
truncating the Fourier series of a CW pulse, as explained in Section III-B-1a. The 
reference signal is a rectangular-envelope CW pulse of duration 20 ms and 
repetition period of 400 ms. The carrier frequency is £ = 500 Hz. Seventeen 
frequency components are used for the transmitted signal, that is, K,,,— 8. The 
resultant bell-shaped CW pulse has a duration of 100 ms. Ringing is barely visible, 
due to the low sidelobes of the Hamming window used for the truncation of the 


Fourier series. The extension of time in the plot of Fig. 21 is equal to one period. 
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Fig. 20. Overall Transfer Function for the Surface Reflection Problem. 
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The acoustic characteristics of the media are the same as for the examples 
in Section III-B-3d and are listed in Table 1. The distance is R$ = 1,000 m. The 
transmitter is located at (Ts, ws 2) = (0, 30, 0) meters, and the receiver at 
(tr, vr, zr) = (1000, 50, 0) meters. 

Figure 22 shows the received pulse when the transmitter is a single point 
source. The distortion due to the multipath nature of the medium is evident. 

Figure 23 shows the received pulse due to a 5x 5 transmitter array. The 


interelement spacing is roughly 4/2 in both the X and Y directions. Note that the 
pulse arrivals occur at the same instants as in Fig. 22, but the relative amplitudes 


are different. 


D. CONCLUSIONS 


The equations that describe a range-independent acoustic communication 
channel were derived by using linear systems theory, a basic engineering tool, as a 
framework. They incorporate aspects of signal processing and acoustic propagation. 
The main results from Chapter Il are the expression for the overall transfer 
function, Eq. (51) or its alternative form, Eq. (60); and the procedure to obtain the 
transformed medium transfer function, the solution to the inhomogeneous 
Helmholtz wave equation, Eq. (90). These results, together with Eq. (40c). were 
used to implement a waveform prediction program, as described in Section III-B. 
We do not claim novelty of results. Equation (60) is a mere extension of a classical 


result, in which the contributions of the point sources that make up an array are 
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Fig. 22. Output Electrical Signal Due to a Single Point Source. 
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simply added, a consequence of the linearity of the wave equation that describes 
small-amplitude acoustic phenomena. The medium transfer function for the layered 
waveguide was derived in order to enable us to implement a working computer 
program. 

The program described in Section III-B is based on a slightly modified form of 
Eq. (60), in which we use the wavenumber y-component k,, as the independent 
variable. The only justification for this change of variable is the easier analysis of 
the medium transfer function. The program works, but is slow. The computation 
of the waveform shown in Fig. 23, for example, took about three minutes per 
frequency component, with an overhead time—independent of the number of 
irequency components—of about two minutes, in an IBM 3033U/4381 2-cpu 
network. Both the speed and accuracy can be improved by using an adaptive 
technique for the truncation of the interval of integration, that is, computation of 
the limits of integration as a function of the behavior of the medium and range 
factors in the integrand. 

Other implementations are possible, for example, the Fast-Field-Program 
described in [Ref. 4:pp. 90-92]. Equation (51). in the form of a two-dimensional 
Fourier transform, could be useful for implementation. 

Concerning the medium transfer function, analytical solutions to the 
range-independent wave equation are only available for a few speed-oí-sound 
profiles. Approximate solutions, like those provided by the WKB method, are 
candidates for implementation. On the other hand, numerical methods may not be 


feasible for our purpose, waveform prediction, due to time constraints. 
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As mentioned in the beginning of Chapter III, the main thrust for the present 
implementation of the range-independent equations was to have a working program. 
The program is working. More than an improvement, it possibly needs a different 
implementation. Nevertheless, the basic important fact is that, independent of the 
implementation, the modular nature of the equations should be fully exploited. 
Solutions to specific problems can be implemented by adding new subprograms, 


with a minimum of code change. Further investigation is recommended. 
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